#!/usr/bin/env python
# encoding: utf-8
"""
pitchfork.py

Created by Andrew M.C. Dawes on 2008-03-21.
Copyright (c) 2008 Andrew M.C. Dawes.
Some rights reserved. License: Creative Commons GNU GPL:
http://creativecommons.org/licenses/GPL/2.0/
"""

import sys
import os
from scipy.integrate import ode
from scipy.io import write_array
import scipy
from scipy import cos, sin, pi
import pylab

save_crossings = True

fig_size = [4,3]
params = {'backend':'ps', 'axes.labelsize':10, 'legend.fontsize':10,'xtick.labelsize':8,'ytick.labelsize':8,'figure.figsize':fig_size}

pylab.rcParams.update(params)

def f(t,y,r,h):
    # return -r*(sin(6*y)+h*sin(2*y))
	return 2*(3+h)*sin(2*y) + 6*r*sin(6*y)
def jac(t,y,r,h):
    # return -r*(-2*h*cos(2*y)+6*cos(6*y))
	return 4*(3+h)*cos(2*y) + 36*r*cos(6*y)

def main():
	y0 = 0.75
	t0 = 0.0
	crossings = []

	for index in range(0,15):
		# h = 10**(-index/5.0)
		h = index*(0.1/15)
		# y0 = index*(pi/15)
		rparam = 1
		r = ode(f,jac).set_integrator('vode').set_f_params(rparam,h).set_jac_params(rparam,h).set_initial_value(y0,t0)
		dt = 0.01
		t1 = 2.0-dt
		tlist = []
		ylist = []
		switch_is_off = True
		while r.successful() and r.t < t1:
			r.integrate(r.t+dt)
			# print r.t, r.y[0]
			if r.y[0] > 0.9 and switch_is_off and save_crossings:
				print "crossing"
				switch_is_off = False
				crossings.append((h,r.t))
			tlist.append(r.t)
			ylist.append(r.y[0])
		
		# print tlist
		pylab.plot(tlist,ylist,"-",lw=0.5,label=("%3.2f" % h))
	crossings_array = scipy.array(crossings)
	pylab.legend(loc="upper right")
	pylab.xlabel('t')
	pylab.ylabel('x')
	pylab.savefig("response_circ.png")
	
	if save_crossings:
		pylab.figure()
		pylab.loglog(crossings_array[:,0],crossings_array[:,1],".")
		pylab.xlabel('h')
		pylab.ylabel('Response time (Arb. Units)')
		# pylab.show()
		pylab.savefig("times_log_circ.png")
	
		pylab.figure()
		pylab.plot(crossings_array[:,0],crossings_array[:,1],".")
		pylab.xlabel('h')
		pylab.ylabel('Response time (Arb. Units)')
		# pylab.show()
		pylab.savefig("times_linear_circ.png")
		write_array("pitchfork_rt_circ.dat",crossings_array)


if __name__ == '__main__':
	main()

